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Abstract 

In  fuel  cell  systems,  a  multitude  of  coupled  physical  and  chemical  processes  take  place  within  the  assembly:  fluid  flow,  diffusion,  charge  and 
heat  transport,  as  well  as  electrochemical  reactions.  For  design  and  optimisation  purposes,  direct  numerical  simulation  of  the  full  three- 
dimensional  (3D)  structure  (using  CFD  tools)  is  often  not  feasible  due  to  the  large  range  of  length  scales  that  are  associated  with  the  various 
physical  and  chemical  phenomena.  However,  since  many  fuel  cell  components  such  as  gas  ducts  or  current  collectors  are  made  of  repetitive 
structures,  volume  averaging  techniques  can  be  employed  to  replace  details  of  the  original  structure  by  their  averaged  counterparts.  In  this 
study,  we  present  simulation  results  for  SOFC  fuel  cells  that  are  based  on  a  two-step  procedure:  first,  for  all  repetitive  structures  detailed  3D 
finite  element  simulations  are  used  to  obtain  effective  parameters  for  the  transport  equations  and  interaction  terms  for  averaged  quantities. 
Bipolar  plates,  for  example,  are  characterised  by  their  porosity  and  permeability  with  respect  to  fluid  flow  and  by  anisotropic  material  tensors 
for  heat  and  charge  transport.  Similarly  one  obtains  effective  values  for  the  Nernst  potential  and  various  kinetic  parameters.  The  complex 
structural  information  is  thereby  cast  into  effective  material  properties.  In  a  second  step,  we  utilise  these  quantities  to  simulate  fuel  cells  in  2D, 
thereby  decreasing  the  computation  time  by  several  orders  of  magnitude.  Depending  on  the  design  and  optimisation  goals,  one  chooses 
appropriate  cuts  perpendicular  or  along  the  stack  axis.  The  resulting  models  provide  current  densities,  temperature  and  species  distributions  as 
well  as  operation  characteristics.  We  tested  our  method  with  the  FEM-based  multiphysics  software  NMSeses,  which  offers  the  flexibility  to 
specify  the  necessary  effective  models.  Results  of  simulation  runs  for  Sulzer  HEXIS-SOFC  stacks  are  presented. 
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1.  Introduction 

Fuel  cells  (FCs)  have  the  ability  to  directly  convert 
hydrogen  or  hydrocarbon  fuels  to  electricity  and  heat  using 
electrochemical  reactions.  Hence,  they  represent  an  efficient 
and  environmentally  friendly  alternative  to  traditional  com¬ 
bustion  processes.  The  commercial  potential  of  fuel  cells  for 
use,  for  example,  in  automotive  applications,  residential  co¬ 
generation,  or  portable  consumer  electronics,  has  been 
recognised  for  many  years.  However,  only  in  recent  years 
major  technical  barriers  have  been  overcome  so  that  FC 
technology  now  has  become  a  major  industrial  trend. 

The  recent  transfer  of  R&D  activities  from  universities  to 
companies  shifted  the  main  focus  from  demonstration  of 
new  materials  and  principles  to  the  development  of  reliable 
and  cost-effective  products  for  mass  market  applications. 
This  trend  also  has  an  impact  on  the  numerical  modelling 
of  FC  devices.  Companies  are  now  in  need  of  effective 
simulation  tools  to  speed  up  their  development  processes. 
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However,  since  fuel  cells  are  complex  devices  with  a  multi¬ 
tude  of  tightly  coupled  physical  and  chemical  phenomena 
taking  place,  it  is  clear  that  there  is  no  simple  simulation  tool 
ready  to  satisfy  all  these  demands.  Rather,  tools  are  needed 
with  a  high  degree  of  flexibility  to  combine  submodels  of 
different  complexity  (incorporating,  for  example,  effects 
occurring  at  different  length  scales)  to  larger  units.  As  an 
example,  consider  the  close  interrelation  of  the  generation  of 
electric  energy  at  the  triple  phase  boundary  with  the  release 
of  heat  due  to  reversible  and  irreversible  processes.  This 
happens  on  a  microscopic  level.  However,  heat  has  to  be 
removed  over  large  distances  to  guarantee  optimal  opera¬ 
tion.  Furthermore,  heat  removal  must  not  ruin  the  electrical 
performance,  e.g.  due  to  large  amounts  of  excess  air  in  the 
cathode  compartment  for  cooling  purposes  or  by  the  need  of 
additional  cooling  fluids. 

Another  feature  of  fuel  cells  systems  is  their  hierarchical 
arrangement  of  many  similar  or  identical  substructures:  a 
FCstack  consists  of  cells,  which  themselves  are  composed  of 
many  similar  flow-field  structures  (meanders,  nipples,  etc.). 
This  enables  one  to  perform  simulations  of  the  smallest  re¬ 
petitive  structures  present  and  to  characterise  their  behaviour 
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through  effective  transport  parameters  such  as  average 
permeabilities  and  average  thermal  conductivities. 

These  detailed  submodels  are  then  included  into  higher- 
level  models  as  geometrically  uniform  materials  with  the 
same  characteristics.  This  approach  is  especially  efficient 
when  a  complex  three-dimensional  (3D)  geometry  can  be 
mapped  onto  a  two-dimensional  (2D)  domain. 

In  this  study,  we  present  such  an  approach  to  simulate 
FCs,  namely  the  numerical  volume  averaging  method 
(NVAM),  which  is  based  on  the  volume  averaging  method 
(VAM)  [1]  and  the  theory  of  macro  transport  equations 
(MTE)  [2].  On  the  one  hand,  this  approach  is  much  more 
efficient  than  the  more  traditional  use  of  CFD  codes  that 
relies  on  the  full  discretisation  of  structural  details  [3]. 
NVAM  enables  the  abstraction  of  a  (part  of  a)  complex 
geometry.  On  the  other  hand,  NVAM  is  rigorously  based  on 
the  governing  balance  and  constitutive  equations  of  con¬ 
tinuum  physics,  i.e.  it  is  not  necessary  to  assess  FC  beha¬ 
viour  by  system  identification  tools  using  a  black  box 
approach.  Therefore,  NVAM  is  able  to  bridge  the  gap 
between  numerical  simulation  on  a  first  principles  level 
and  a  pure  system  dynamics  modelling  approach. 

2.  Averaged  transport  equations 

To  illustrate  our  approach  we  consider  as  an  example  the 
mass  transport  in  porous  media.  The  term  porous  is  used  in  a 
wide  sense,  here  we  consider  the  gas  diffusion  in  an  SOFC 
electrode  as  well  as  the  flow  field  in  millimetre  range 
channels  as  porous  flow.  The  method  then  can  be  generalised 
to  additional  transport  phenomena  of  interest  in  FC  model¬ 
ling,  including  heat  and  charge  transport,  etc. 

2.7.  The  volume  averaging  method 

The  volume  averaging  method  (VAM)  was  developed  to 
describe  transport  phenomena  in  ordered  or  disordered 
porous  media  [1].  The  existence  of  at  least  two  different 
length  scales  is  characteristic  for  these  media:  the  typical 
diameter  of  the  pores  and  the  device  length,  as  shown  in 
Fig.  1.  The  actual  transport  process  takes  place  in  the  void 

L 

< - > 


Fig.  1.  Length  scales  in  a  porous  medium:  L,  device  length  scale,  /p  pore 
diameter.  The  radius  r  of  the  averaging  volume  has  to  satisfy  the  relation 
/p  r  L. 


space  defined  by  the  pores.  The  governing  equations  are  due 
to  Navier-Stokes,  i.e. 

V(pv  0  v)  =  —  VP  +  V  •  t 

Xij  =  fi(diVj  +  djVj)  - 1  fi(v  ■  v)Sjj,  (1) 

with  the  pressure  P,  the  fluid  density  p,  the  viscosity  p,  and 
the  stress  tensor  t.  These  quantities  are  defined  only  within 
the  pores.  Eq.  (1)  are  virtually  intractable  numerically  for 
many  technically  important  structures  due  to  the  complex 
geometry  and  the  immense  effort  to  specify  a  computational 
grid  and  the  boundary  conditions  for  a  piece  of,  e.g.  a  porous 
ceramic  material. 

The  solution  to  this  problem  is  the  description  of  physi¬ 
cally  relevant  fields  only  on  an  average  basis.  The  approach 
given  in  [1]  deduces  effective  transport  equations  for  sui¬ 
tably  defined  volume  averages  of  the  true  physical  quan¬ 
tities:  the  averaged  velocity  field  is  defined  by 

Vav  =  T  J J  J  v(x)dV.  (2) 

where  Q  denotes  an  averaging  volume  that  is  large  with 
respect  to  the  pore  length  scale,  but  small  compared  to  the 
device  length.  Vq  is  the  volume  of  Q.  After  a  lengthy 
calculation  one  finds  a  modified  transport  equation  for  the 
averaged  velocity  field  together  with  an  updated  expression 
for  the  stress  tensor  Teff 

-VP  +  V  •  ?eff  -  7 Vav  =  0,  (3) 

k 

where  k  denotes  the  permeability  of  the  porous  material. 
Usually  the  term  involving  the  stress  tensor  is  negligible 
compared  to  the  last  one,  leading  after  rearrangement  to 

7  vav  =  -  VP  (4) 

k 

the  well-known  Darcy  equation.  This  equation  applies  to  the 
whole  domain  of  the  device,  i.e.  there  is  no  longer  a  need  to 
formulate  boundary  conditions  within  the  pores  or  to  dis- 
cretise  the  pore  geometry  as  the  averaged  velocity  field  is 
continuous  across  the  whole  domain.  The  quantity  k,  which 
is  a  tensor  of  rank  two  in  general,  now  contains  the  informa¬ 
tion  about  the  pore  geometry  and  expresses  an  anisotropy 
due  to,  e.g.  an  alignment  of  the  pores  along  a  characteristic 
direction  of  the  material. 

The  challenge  in  applying  VAM  to  transport  problems  is 
the  solution  of  the  so-called  closure  equations  [1],  which 
implicitly  define  the  effective  parameters  characterising  the 
porous  material,  (p/k  in  Eq.  (4)).  The  closure  is  a  partial 
differential  equation  with  at  least  the  same  complexity  as  the 
original  transport  equation.  It  has  to  be  solved  on  a  repre¬ 
sentative  region  of  the  porous  medium  using  suitable  bound¬ 
ary  conditions.  Sometimes  it  is  possible  to  deduce  certain 
properties  of  the  effective  parameters  using  analytical  meth¬ 
ods,  but  this  rarely  leads  to  actual  figures  for  the  transport 
parameters. 
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Fig.  2.  Periodic  porous  structure  and  a  possible  elementary  cell. 


2.2.  The  macro  transport  equation  method 

In  [2],  a  slightly  different  approach  is  outlined  to  deal  with 
transport  in  porous  media:  instead  of  starting  with  the 
definition  of  volume  averages  for  the  relevant  quantities 
and  deducing  their  transport  equations,  so-called  macro 
transport  equations  (MTE)  are  set-up.  These  equations  allow 
for  the  calculation  of  the  physical  quantities  at  suitably 
chosen  reference  points.  The  method  is  ideal  for  ordered 
materials,  but  also  applies  with  certain  restrictions  to  dis¬ 
ordered  materials. 

One  of  the  critical  issues  with  this  method  is  the  deter¬ 
mination  of  the  effective  transport  parameters.  Again  one 
has  to  solve  specific  abstract  PDEs  on  elementary  cells  on  a 
periodic  domain,  cf.  Fig.  2. 

2.3.  The  numerical  volume  averaging  method 

To  simplify  the  application  of  either  VAM  or  MTE  one 
has  to  avoid  the  solution  of  the  closure  PDEs.  Instead  of 
doing  abstract  mathematics,  we  perform  virtual  experiments 
on  a  repetitive  region  of  the  porous  medium  using  a  com¬ 
puter  simulation  code  to  obtain  the  effective  parameters,  as 
shown  in  Fig.  3.  With  suitable  boundary  conditions  on  this 
domain  we  are  able  to  “measure”  the  tensor  components  of 
the  permeability  tensor  k ,  i.e.  kxx  and  kyy.  This  is,  of  course, 
reminiscent  to  the  work  of  an  experimentalist,  who  does 
exactly  the  same  thing  in  reality  to  confirm  Darcy’s  law  for  a 
given  material.  He  sets  up  a  pressure  gradient  in  a  chosen 
direction,  measures  the  corresponding  flow  rates  in  ah 
possible  directions,  and  extracts  the  material  properties  in 
tensor  form. 

Although  NVAM  still  is  based  on  the  solution  of  PDEs  on 
complex  domains,  it  has  a  number  of  advantages  over  its 
parent  methods: 

•  There  is  no  need  to  implement  specific  closure  PDEs  as  in 
VAM  to  get  effective  parameters.  One  uses  the  given,  well 
known  transport  equations  with  their  clear  physical  inter¬ 
pretation. 


Fig.  3.  NVAM  for  a  complex  pore  structure  showing  anisotropic 
behaviour:  calculation  of  kxx  (top)  and  of  kyy  (bottom).  Solid  lines  indicate 
the  driving  boundary  conditions  for  the  flow,  while  on  the  remaining 
boundaries  the  no-flux  condition  applies. 

•  There  exists  a  multitude  of  commercial  and  non-com¬ 
mercial  numerical  simulation  tools  that  already  perform 
the  necessary  calculations. 

•  With  the  analogy  to  experimental  work  and  the  physical 
interpretation  of  the  necessary  simulations,  it  is  possible 
to  apply  NVAM  to  virtually  any  transport  phenomenon. 

•  NVAM  also  works  for  complex  material  behaviour  found 
in  many  engineering  applications.  Consider,  e.g.  the  case, 
where  mass  transport  occurs  in  non-isothermal  states. 
While  VAM  in  its  basic  formulation  is  not  applicable 
because  isothermal  behaviour  is  assumed,  NVAM  still 
works:  one  performs  the  virtual  experiment  simply  by 
including  the  given  thermal  gradient.” 

•  Moreover,  it  is  possible  to  extend  the  approach  to  more 
complex  situations,  where  transport  processes  take  place 
in  conjunction  with  chemical  reactions.  While  the  gen¬ 
eralisation  of  VAM  or  MTE  is  complex,  there  is  a 
straightforward  generalisation  in  terms  of  NVAM. 

•  VAM  or  MTE  are  used  to  rigorously  deduce  the  form  of 
the  effective  transport  equations,  while  NVAM  provides  a 
framework  for  efficient  calculations  for  engineering  tasks. 

2.4.  Numerical  algorithms 

In  principle  any  numerical  discretisation  scheme  could  be 
adopted  for  the  actual  solution  of  mass  transport  or  diffusion 
equations,  etc.  The  most  prominent  example  is  the  finite 
volume  method  in  conjunction  with  an  explicit  time  inte¬ 
gration  scheme.  Another  choice  is  the  finite  element  method 
(FEM),  suitable  especially  for  stationary  states.  We  have 


1  The  macro  transport  equation  method  is,  in  fact,  much  more  general  2  However,  it  is  no  longer  guaranteed  that  the  form  of  the  effective 

and  allows  also  for  the  simplification  of  transport  equations  if  the  transport  equation  still  remains  valid.  This  could  be  proven  only  by  regress 

simulation  domain  contains  a  direction  with  reduced  complexity,  e.g.  a  slit.  to  the  traditional  approach,  i.e.  VAM. 


M.  Roos  et  al.  /  Journal  of  Power  Sources  118  (2003)  86-95 


89 


used  this  method  because  it  provides  the  highest  flexibility  to 
simulate  complex  material  behaviour  and  boundary  condi¬ 
tions.  Through  the  volume  averaging  process,  all  the  infor¬ 
mation  about  the  geometrical  structure  is  cast  into  the 
anisotropy  of  the  effective  parameters.  For  layered  media, 
for  example,  different  values  of  the  thermal  conductivity  in 
or  transversal  to  the  layer  result  [2].  Anisotropic  material 
behaviour  can  easily  be  adopted  in  finite  element  formula¬ 
tions.  The  same  is  true  for  complex  interaction  models, 
which  arise,  e.g.  by  averaging  the  electrochemical  kinetics 
over  electrode  and  electrolyte  domains.  As  a  result,  the 
phenomenological  kinetics  depend  on  the  current  density 
due  to  mass  transport  limitations  through  the  porous  elec¬ 
trode. 

2.5.  Application  to  fuel  cell  structures 

We  applied  NVAM  to  fuel  cell  simulations.  This  is  of 
great  interest  since  a  typical  fuel  cell  device  contains  at  least 
four  different  length  scales  rendering  numerical  simulations 
especially  cumbersome.  These  are  (with  an  example  of  the 
dominant  effect): 

_ rj  _ o 

•  Nanopore  scale  (from  10  to  10  m):  chemical  reac¬ 

tions  occurring  at  the  triple  phase  boundary. 

•  Porous  material  scale  (from  10  to  10  m):  diffusion 

transport  in  ceramic  electrodes. 

•  Flow  channel  scale  (from  10-3  to  10-4  m):  free  fluid  flow 
together  with  diffusion. 

•  Device  scale  (from  10  to  10  m):  cells  and  stacks 

operation,  including  the  balance  of  plant. 

Of  course,  it  is  impossible  with  desktop  computers  to 
cover  this  range  of  length  scales  within  a  single  numerical 
scheme.  The  macro  transport  equations  approach  can  serve 
as  a  generic  method  to  tackle  this  issue.  One  starts  from  the 
physical  description  level  and  develops  a  cascade  of  macro 
transport  equations:  each  level’s  transport  equations  are 
parameterised  by  models  of  the  next  lower  level.  Put  into 
action,  this  scheme  works  out  as  follows  for  SOFC  simula¬ 
tion: 

•  A  3D  model  of  the  micropores  at  the  electrode-electrolyte 
interface  provides  effective  kinetic  models  for  the  elec¬ 
trochemical  processes  at  the  next  higher  level. 

•  A  3D  model  of  a  representative  part  of  a  MIC/electrode/ 
electrolyte/MIC  structure  is  used  to  obtain  effective  trans¬ 
port  parameters  for  fluid  transport  at  the  flow  channel 
scale  in  the  MIC. 

•  A  2D  model  of  the  cell  (or  the  stack)  is  set-up.  To  describe 
fluid  flow,  the  MIC  is  modelled  as  a  macroporous  structure. 

•  Stacking  of  several  cell  models  allows  for  the  simulation 
of  whole  stacks. 

It  is  possible  to  combine  this  approach  with  experimental 
data  at  any  level.  If  the  lowest  level  is  omitted,  missing 
information  can  be  taken  from  experimental  work,  e.g. 
impedance  spectroscopy  to  access  the  basic  parameters  of 


the  electrochemical  kinetics.  This  is  often  necessary,  as  the 
elementary  chemical  reactions  are  rarely  known.  At  any  level, 
the  experimental  data  can  be  used  to  validate  NVAM  results. 

2.(5.  Lateral  models 


The  typical  geometry  of  flow  fields  in  PEFC  cells  is  truly 
three-dimensional.  With  approximate  2D  models  of  these 
structures,  some  relevant  interactions  or  geometric  influ¬ 
ences  have  to  be  discarded.  NVAM  provides  a  novel 
approach  to  this  issue:  a  lateral  2D  cell  model  is  specified 
by  adopting  suitable  macro  transport  equations  for  an  aver¬ 
aged  flow  field.  The  omitted  third  dimension  (z-direction  in 
Fig.  4)  manifests  itself  in  various  ways:  the  viscous  friction 
with  the  upper  and  lower  walls  is  absorbed  into  the  perme¬ 
ability  tensor.  Averaging  according  to  NVAM  provides  the 
anisotropic  tensor  for  the  permeability.  In-plane  heat  transfer 
is  also  described  by  anisotropic  material  tensors.  Heat  and/or 
charge  transport  to  neighbouring  cells  in  the  transverse  (z- 
direction)  is  accounted  for  by  setting  up  suitable  interaction 
models,  e.g.  a  temperature  dependent  heat  source  to  reflect 
Newton-type  cooling. 

The  electrochemical  reactions  (taking  place  at  an  inter¬ 
face  between  layers)  are  described  in  a  lateral  model  by  an 
effective  2D  volume  reaction.  A  detailed  3D  submodel 
similar  to,  e.g.  the  one  presented  in  [4]  for  PEFC  processes, 
provides  effective  parameters.  Geometric  effects,  properties 
of  the  ME  A,  transversal  diffusion  in  the  gas  channels,  etc.  is 
accounted  for  by  complex  kinetic  equations,  to  be  parame¬ 
terised  with  help  of  the  3D  submodel. 

Due  to  the  reduction  of  complexity  by  using  a  2D  lateral 
model,  localised  effects,  such  as  the  formation  of  hot  spots, 
cannot  be  studied  in  this  way.  Only  phenomena  on  the  length 
scale  of  the  averaging  process  used  in  NVAM  can  be 
addressed. 
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Fig.  4.  The  3D  flow-field  structure  and  its  representation  in  a  2D  lateral 
model  with  anisotropic  material  behaviour. 
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2.7.  Concluding  remarks 

The  description  of  the  NVAM  closes  with  some  remarks 

about  this  method: 

•  The  application  of  the  method  to  actual  devices  is  straight¬ 
forward,  as  it  does  not  depend  on  complex  mathematical 
derivations  or  the  need  for  an  implementation  of  unfa¬ 
miliar  abstract  PDEs. 

•  NVAM  models  can  be  tailored  to  fit  to  engineering  tasks. 
The  researcher  decides  on  the  degree  of  complexity 
necessary  to  tackle  the  actual  problem,  while  keeping 
the  computational  costs  low. 

•  It  is  possible  to  introduce  experimental  data  at  any  levels 
of  the  modelling  hierarchy,  thereby  replacing  the  lower 
levels  or  validate  the  numerical  model. 

•  NVAM  interpolates  from  the  pure  scientific  approach 
with  its  ab  initio  transport  process  formulation  up  to 
the  pure  engineering  task  of  black  box  modelling. 

•  Application  of  NAVM  exhibits  a  strong  analogy  to  experi¬ 
mental  work.  In  fact,  this  modelling  approach  is  best 
described  as  performing  virtual  experiments.  It  is  hoped 
that  it  will  be  accepted  more  easily  by  development 
engineers  than  more  abstract  numerical  schemes. 


3.  Application  to  SOFC  modelling 

We  tested  the  NVAM  by  simulating  solid  oxide  fuel  cells 

a 

and  stacks  with  the  HEXIS  design  [5-7].  This  system 
exhibits  a  high  integration  level,  as  heat  exchange,  fuel 
reforming,  electrochemical  reactions  and  after  burning  take 
place  within  a  single  housing.  The  design  of  the  stack  asks 
for  careful  optimisation  to  meet  the  (often)  conflicting  goals. 
This  can  be  done  efficiently  only  with  help  of  detailed 
numerical  simulations  because  many  processes  strongly 
interact.  Due  to  the  geometrical  complexity  of  the  real 
structure,  only  a  method  like  NVAM  can  support  engineer¬ 
ing  tasks  at  affordable  costs  [7]. 

To  set-up  a  single  cell  model,  we  started  from  the  second- 
lowest  level  (macroporous  scale)  because  there  is  no  detailed 
information  on  the  nano  scale  available.  On  any  stage,  the 
relevant  set  of  PDEs  describing  the  transport  phenomena 
have  to  be  solved.  The  following  set  of  independent  fields 
has  been  used:  temperature  T,  velocity  field  v,  pressure  P, 
molar  fractions  for  each  species  out  of  the  set 
a  e  {H2,  H20, 02,  N2,  CO,  CO2,  CH4},  and  the  electric 
potential  0. 

3.1.  Transport  equations  for  SOFC  systems 


flow-field  regions  of  a  3D  model  of  the  MIC,  or  to  the  piping. 
The  ideal  gas  law  links  the  mass  density  to  the  species  molar 
fractions,  pressure  and  the  local  temperature  by 


with  the  universal  gas  constant  R ,  the  molar  masses  Ma  and 
the  molar  species  fractions  va.  Of  course,  more  complex 
constitutive  laws  for  the  fluid  are  feasible,  but  the  ideal  gas 
law  describes  the  behaviour  of  high  temperature  gases  very 
well.  In  porous  media,  e.g.  within  the  ceramic  electrodes  or 
in  2D  representations  of  3D  flow  fields,  we  adopt  Darcy’s 
law  (Eq.  (3)). 

In  any  case,  one  has  to  combine  Eq.  (1)  or  Eq.  (2)  with  the 
continuity  equation,  i.e.  the  mass  balance  reading  for  sta¬ 
tionary  states 


V  •  (pv)  =  0.  (6) 

The  diffusion  equation  for  species  transport  in  void  and 
porous  materials  reads 

v  •  (pav)  =  -v  7a  +  na,  (7) 

where  the  mass  density  for  species  a  is  given  by  pa  =  vaMac 

and  the  molar  concentration  equals  c  =  p  (^avaMa)  .  77a 
denotes  the  production  rate  of  species  a  (based  on  the  mass 
density).  This  rate  is  the  sum  of  all  active  chemical  reactions, 
i.e. 

nri  =  Tfp  +  77™  (8) 

The  first  term  is  the  contribution  of  all  homogenous  and 
quasi  homogenous  reactions,  e.g.  the  reforming  and  shift 
reactions  in  the  Ni-cermet  anode  of  a  SOFC  system  and 
reads  77^H  =  MaJ]nvanrw,  where  vaw  is  the  stoichiometric 
coefficient  for  species  a  in  reaction  n  and  rn  the  molar 
reaction  rate  density  for  reaction  n.  This  last  quantity  has 
to  be  supplied  according  to  the  catalytic  properties  of  the 
given  material.  Of  course,  the  rate  densities  usually  depend 
on  the  local  temperature,  species  concentrations,  etc.  The 
second  term  on  the  RHS  of  Eq.  (8)  represents  the  contribu¬ 
tions  of  the  electrochemical  reactions  at  the  triple  phase 
boundary  (TPB).  Their  form  will  be  discussed  later. 

The  species  diffusion  current  density  (measured  with 
respect  to  the  centre  of  mass  system),  finally,  reads 


c2MaMp 

P 


Oafi  V . xa , 


(9a) 


where  Stefan-Maxwell’s  law  implicitly  defines  the  effective 
diffusion  constants  Dap  by 


The  governing  law  for  fluid  flow  in  void  regions  are 
the  Navier-Stokes  equations,  Eq.  (1),  which  apply  to  the 


(9b) 


3  The  acronym  HEXIS  corresponds  to  heat  exchanger  integrated  stack,  a 
SOFC  concept  developed  by  Sulzer  AG,  Switzerland. 


The  latter  expression  depends  on  the  binary  diffusion  coef¬ 
ficients  Dap,  which  are  functions  of  temperature  and 
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pressure.  In  microporous  regions,  the  binary  diffusion  coef¬ 
ficients  are  further  scaled  to  account  for  the  effect  of  pore 
diffusion  (Bruggemann  correction).  The  effective  value 
reads  then 

Da[J:p  =  Gh5Dap,  (10) 

where  s  denotes  the  (volumetric)  porosity.  The  heat  transport 
equation  applies  to  all  simulation  domains  (with  vanishing 
convective  contributions  in  solid  materials,  of  course) 


3.2.  Implementation  of  the  electrochemistry 

In  addition  to  (quasi)  homogenous  reactions  in  the  porous 
electrodes  or  void  spaces,  one  has  to  describe  the  electro¬ 
chemical  reactions.  SOFC  systems  are  able  to  electroche- 
mically  oxidise  hydrogen  as  well  as  carbon  monoxide 
directly.  The  overall  reactions  reads 

2H2  +  02  ^  2H20,  (16a) 

and 


v  •  I  I  =  -v  -jth  -  v  •  yfha]a  +  nth,  (i  i) 

\  a  /  a 

where  ha  is  the  (mass  density  based)  enthalpy  of  species  a, 
approximated  for  ideal  gases  by  ha  =  cPAT.  The  heat  pro¬ 
duction  rate  77th  is  given  by 

nt  h  =  Y.nZ  +  E77™  +  +  nlh£XL ,  (12) 

n  s 

where  the  first  term  on  the  RHS  corresponds  to  the  heat 
release  of  the  quasi-homogenous  reactions.  We  set  77^  = 
—A Hnrn  with  A Hn  denoting  the  reaction  enthalpy  for  reac¬ 
tion  n  and  rn  the  corresponding  rate.  The  second  term 
corresponds  to  the  electrochemical  reactions,  see  Section 
3.2. 

1  2 

The  third  term,  7Tth,ei  =  \jq\  ,  is  due  to  Joule’s  heating 
with  the  material’s  electric  conductivity  o  and  the  electric 
charge-current  density  jq.  The  term  77th,ext  accounts  for  an 
external  heat  source  or  sink  due  to  cooling  mechanisms  or — 
for  2D  lateral  models — accounts  for  heat-currents  directed 
into  the  omitted  third  direction,  cf.  remarks  at  the  end  of 
Section  2.7. 

Finally,  the  heat-current  density  reads  according  to  Four¬ 
ier’s  law  of  heat  conduction 

Jth  =  -2VT.  (13) 

The  thermal  conductivity  2  might  depend  on  local  species 
concentrations,  temperature  and  is  in  general  a  tensor  quan¬ 
tity  as  explained  in  Section  2. 

The  charge  transport  equation,  defined  only  in  electrically 
conducting  materials  reads 

-V  -l,  +  E^JPB  =  0,  (14) 


2CO  +  02  =>  2C02,  (16b) 

The  actual  reactions  at  the  anode  triple  phase  boundary  are 
given  by 

2H2  +  202“  =■  2H20  +  4e“  ( 1 7a) 

and 


2CO  +  202“  2C02  +  4e“  (17b) 

while  the  process  at  the  cathode  TPB  reads 

02  +  4e_  =>  202- .  (18) 


These  reactions  contribute  to  the  source  terms  in  the  charge, 
species  and  heat  transport  equations.  The  source  term  for 
charge  transport  reads 

£jPB  =  A^HV(<S(/)-|V/|),  (19) 

where  the  scalar  function/(x)  =  0  defines  the  location  of  the 
TPB  interface  in  space.  AWS  denotes  the  (irreversible) 
potential  jump  across  the  TPB  for  reaction  s  and  is  given  by 


AYS  =  A<F° 


(20) 


where  AT®  denotes  Nemst’s  potential  at  thermodynamic 
equilibrium  for  reaction  s  (cf.  Eq.  (24)),  gs  is  the  double  layer 
surface  conductance,  which  is  used  to  model  the  activation 
overpotential,  cf.  Eq.  (25),  and  jqyS  finally  corresponds  to  the 
current  density  across  the  TPB  due  to  reaction  s. 

The  additional  contributions  to  the  source  terms  for  heat 
and  species  transport  can  be  expressed  with  help  of  the 
reaction  rate  per  unit  area  at  the  TPB 


.TPB  _  (jq,s  *  W)S(f) 


(21) 


where  Es  represents  the  (singular)  charge  double  layer,  cf. 
Section  3.2.  The  charge-current  density  jq  is  given  by  the 
microscopic  version  of  Ohm’s  law 

jq  =  -<7V4>,  (15) 


where  ns  corresponds  to  the  number  of  electrons  transferred 
in  reaction  s  and  F  is  Faraday’s  constant.  The  source  term  in 
Eq.  (12)  due  to  electrochemistry  reads  then 


(22) 


where  a  denotes  the  electric  conductivity  in  MIC  and 
electrode  regions  and  the  ionic  conductivity  in  the  electro¬ 
lyte  region.  It  is  well  known  that  the  ionic  conduction  in 
solid  oxide  electrolytes  can  be  approximated  by  Ohm’s  law 
[10]  (in  contrast  to  liquid  or  polymer  electrolytes). 


where  va?5  is  again  the  stoichiometric  coefficient.  The  con¬ 
tribution  to  the  heat  balance  is  given  by 


rjTPB  _ 
11  th  — 


-r^AA,rJPB  +  E^-1 1«  -JqA2m\9f\,  (23) 
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where  ASs  denotes  the  entropy  release  in  reaction  s  (rever¬ 
sible  release  of  heat)  and  the  second  term  corresponds  to 
losses  due  to  irreversible  processes,  i.e.  the  activation  losses. 

The  expression  for  Nernst’s  potential  for  reaction  s  at 
thermodynamical  equilibrium  reads 


A  ¥?  = 


A Hs  -  TSSs 
nsF 


RT  ,  - 

- >  v 

n*F  ' 


a  ,s 


ln(xa) 


a 


(24) 


and  Block2  =  {species  molar  fractions,  temperature,  elec¬ 
tric  potential}. 

The  overall  loop  converges  not  exceeding  15  Newton- 
Raphson-iterations,  taking  typically  <2  min  computer  time 
per  SOFC  operation  point  on  an  Athlon  XP2000  (1.6  GHz) 
for  the  single  cell  model  of  Section  3.5.  Therefore,  a  full 
performance  plot  is  calculated  in  <2  h.  The  3D  models  of 
Section  3.4  generally  require  more  computer  resources. 


where  the  first  term  corresponds  to  Nernst’s  potential  at 
normal  conditions  and  the  second  contribution  accounts  for 
the  change  of  the  activity  in  the  actual  state.  If  there  is  more 
than  one  electrochemical  reaction  present  (e.g.  the  oxida¬ 
tion  of  CO  in  the  solid  oxide  fuel  cell)  the  actual  potential 
jump  at  the  TPB  is  determined  by  the  competing  reactions. 
Even  for  zero  net  current  density  across  the  electrolyte, 
there  might  be  non-zero  rates  for  the  individual  reactions.  It 
is  assumed  here,  that  they  interact  only  via  their  contribu¬ 
tion  to  the  total  potential  jump,  but  can  be  treated  indepen¬ 
dently  otherwise. 

The  activation  overpotential  at  cathode  or  anode  is  impli¬ 
citly  defined  through  the  surface  conductance  [8]  and  reads 


772  f 


IF 

RT 


(Eact,s)/(RT) 


(25) 


Here  ks  denotes  a  phenomenological  pre-exponential  factor 
for  reaction  s,  Eact  s  its  the  activation  energy,  Pref  a  reference 
pressure,  and  ms  an  exponent  describing  the  activation 
kinetics.  This  corresponds  to  a  linearisation  of  the  Butler- 
Volmer  equation  and  is  well  justified  for  SOFC  conditions 
[8]. 


3.3.  Solution  method 


The  set  of  partial  differential  equations,  i.e.,  Eqs.  (1)  or 
(4),  (6),  (7),  (9a),  (9b),  (11),  (13)-(15)  have  been  imple¬ 
mented  in  the  FE-based  code  NMSeses  [9].  Most  material 
and  process  properties  can  be  accessed  on  the  user  level,  i.e. 
material  temperature  dependencies,  kinetic  laws,  etc.  can  be 
freely  defined. 

Due  to  the  strong  non-linear  dependencies  in  the  expres¬ 
sions  for  the  reaction  rates  and  some  of  the  material  proper¬ 
ties,  the  FE-tool  allows  for  the  specification  of  partial 
derivatives  (with  respect  to  the  independent  fields)  to  speed 
up  the  convergence  rate. 

Grouping  the  independent  fields  into  two  or  more  solution 
field  blocks  can  further  optimise  the  convergence  behaviour. 
The  fields  within  a  given  block  are  solved  by  the  Newton- 
Raphson  algorithm,  which  converges  quadratically  if  the 
initial  starting  solution  is  carefully  chosen  and  the  partial 
derivatives  w.r.t.  all  block  fields  are  defined. 

The  algorithm  implemented  in  NMSeses  solves  the  given 
blocks  consecutively  until  overall  convergence  is  obtained. 
This  is  usually  fast,  if  an  optimal  block  structure  has  been 
defined.  Best  convergence  properties  have  been  found 
with  a  two-blocks  choice:  Blockl  =  {velocity,  pressure} 


3.4.  3D  model  at  macroporous  level 
3.4.1.  Model  specification 

On  the  macroporous  level,  we  used  a  3D  model  for  the 
domain  shown  in  Fig.  5.  It  reflects  the  true  geometry  of  the 
flow  regions  inside  the  MIC  described  by  Eq.  (1).  For 
cathode  and  anode  electrodes,  the  porous  flow  equation 
(Eq.  (4))  is  used.  The  material  properties,  i.e.  porosity 
and  permeability  were  taken  from  the  literature  [10],  or 
fitted  to  experimental  results  if  available,  or  estimated 
otherwise. 

An  important  issue  in  applying  NVAM  is  the  specification 
of  boundary  conditions.  This  is  straightforward  for  heat 
conduction  or  laminar  flow:  the  symmetries  of  the  device 
dictate  the  boundary  conditions  of  the  domain.  We  have  used 
zero-flow  conditions  across  the  mirror  planes,  while  planes 
perpendicular  to  the  flow  direction  were  modelled  isother¬ 
mal  (or  isobaric)  with  a  given  drop  of  temperature  (or 
pressure)  across  the  system. 

Setting  up  the  boundary  conditions  becomes  more 
involved  for  systems  where  dispersion,  i.e.  the  interaction 
of  diffusive  and  convective  transport  (of  heat  or  species) 
plays  a  major  role.  In  this  case,  the  boundary  conditions  are 
less  obvious  because  there  are  fewer  symmetries  to  use,  as 
the  direction  of  the  convective  current  breaks  one  symmetry. 


Fig.  5.  Side  (top)  and  horizontal  (bottom)  view  of  the  metallic 
interconnect  (MIC)  structure.  The  repetitive  element  used  for  numerical 
simulation  is  indicated  by  the  dashed  lines. 
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This  can  be  tackled  by  adopting  generalised  periodic  bound¬ 
ary  conditions.  While  the  heat-current  density  (or  the  velo¬ 
city  flow)  is  periodic  across  the  repetitive  domain  in  flow 
direction,  the  corresponding  potential  field,  i.e.  temperature 
(or  pressure)  consists  of  a  periodic  part  plus  a  superimposed 
linear  ramp  [2].  In  the  parameter  range  for  the  actual 
simulations,  we  were  able  to  estimate  the  transport  para¬ 
meters  for  heat  and  fluid  flow  independently,  neglecting 
dispersion. 

To  date  no  simulations  with  electrochemical  reactions 
have  been  performed  with  the  3D  model  shown  in  Fig.  5. 
While  such  calculations  are  of  great  interest  to  predict  the 
net  efficiency  of  the  energy  conversion  and  their  dependence 
on  geometry,  it  is  mandatory  to  have  independent  measure¬ 
ments  of  the  local  properties  of  the  electrochemical  inter¬ 
actions.  Since  so  far  only  IV-curves  for  full  cells  are 
available,  no  direct  comparison  of  NVAM  and  experiment 
is  possible. 

3.4.2.  Results  for  the  3D  model 

We  calculated  the  permeability  for  the  flow  fields  in  the 
MIC  structure  for  the  anode  and  cathode,  as  well  as  the  heat 
and  charge  conductivities.  Unfortunately,  there  are  no  inde¬ 
pendent  measurements  for  these  transport  parameters  avail¬ 
able.  Nevertheless,  the  numerical  simulations  still  provide 
important  insight:  one  is  able  to  estimate  trends  if  the 
geometrical  structure  is  changed  under  otherwise  identical 
conditions.  These  trends  are  often  much  more  reliable  in 
comparison  to  the  predictions  of  absolute  values. 

3.5.  Effective  2D  cell  model 

The  transport  parameters  from  the  3D  submodel  (see 
Section  3.4)  were  readily  included  into  the  rotational  sym¬ 
metric  2D  model,  see  lower  part  in  Fig.  6.  The  reaction 
properties  have  to  be  taken  from  experimental  IV-curves 
with  help  of  the  Ansatz  functions  for  SOFC  electrochem¬ 
istry,  cf.  Eq.  (22). 

3.5.1.  Model  specification 

We  have  chosen  a  rotational  symmetric  domain  to 
approximate  the  HEXIS  cell.  The  process  of  NVAM 


Fig.  6.  The  top  view  shows  a  part  of  a  HEXIS-SOFC  MIC,  while  the 
lower  sketch  depicts  the  corresponding  2D  rotational  symmetric  FE  model. 


Fig.  7.  Molar  fractions  for  the  species  H2,  H20,  CO  and  C02  along  the 
radial  coordinate  of  the  2D  cell  model  shown  in  Fig.  6. 


(with  an  averaging  volume  large  compared  to  the  exten¬ 
sion  of  the  nipples)  renders  the  MIC  flow-field  rota¬ 
tional  symmetric.  This  is  not  true,  however,  for  the  after 
burning  and  the  air  inlet  zones  along  the  outer  circum¬ 
ference.  At  the  present  stage  we  have  to  accept  this 
approximation. 

Fig.  7  shows  the  molar  fractions  of  hydrogen,  water, 
carbon  monoxide,  and  carbon  dioxide  along  the  radial  co¬ 
ordinate.  The  rather  steep  decrease  of  the  gas  components 
near  60  mm  is  due  to  the  onset  of  the  after-burning  reaction 
and  the  diffusive  mixing  with  the  species  of  the  air  gas 
stream.  Interpretation  of  these  model  properties  under 
varying  operation  conditions,  such  as  fuel  mass  flow,  stack 
temperature,  etc.,  is  the  key  to  improve  the  understanding 
of  the  internal  states  of  SOFC  cells.  Since  direct  measure¬ 
ments  of  the  gas  composition  inside  the  cell  is  very 
difficult,  these  results  can  be  validated  only  indirectly, 
e.g.  by  comparing  IV-curves  or  performance  maps,  cf. 
Section  3.5.2. 

3.5.2.  IV- Curves  and  performance  maps 

After  the  model  calibration  using  a  data  set  for  a  single 
value  of  the  fuel  mass  flux,  the  model  can  be  used  to  simulate 
the  performance  at  different  mass  fluxes  and  external  loads. 
In  Fig.  8,  a  typical  set  of  IV-curves  is  shown,  for  varying  fuel 
mass  fluxes. 

This  two-parameter  set  of  results  can  be  cast  into  so-called 
performance  maps,  Fig.  9.  It  shows  lines  of  constant  effi¬ 
ciency  (solid  lines),  electrical  power  (dashed)  and  fuel 
utilisation  (dashed-dotted).  This  information  is  plotted  in 
a  fuel  mass  flow  versus  voltage  graph.  The  two  graphs  in 
Fig.  9  correspond  to  a  thick  electrolyte  (top  graph)  and  thin 
one  (bottom),  respectively.  A  comparison  clearly  demon¬ 
strates  improved  performance  of  the  cell  with  the  thinner 
electrolyte.  Many  more  variations  of  the  modelling  condi¬ 
tions  are  feasible.  The  extraction  of  performance  trends  with 
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Fig.  8.  IV-curves  for  different  methane  mass  fluxes,  1  g/h  (solid),  3  g/h 
(dashed)  and  5  g/h  (dashed-dotted). 


respect  to  maximum  power  or  efficiency  help  to  improve  a 
redesign  of  SOFC  stacks. 

3.6.  System  models 

A  further  step  in  the  modelling  hierarchy  addresses  the 
system  level.  In  order  to  maintain  a  nearly  constant  tem¬ 
perature  throughout  the  stack  at  different  operating  condi¬ 
tions,  the  thermal  insulation  around  the  stack  has  to  be  able 
to  dynamically  adjust  to  the  amount  of  heat  released  in  the 
after-burning  zone.  This  can  be  achieved  through  control  of 
the  air  flow  rate  from  the  outside  to  the  inside  of  the 
insulation  and  through  effective  heat-exchange  between 
the  cold  air  and  the  hot  exhaust  gases. 

Fig.  10  shows  the  smallest  repetitive  structure  of  an 
arrangement  of  two  insulation  walls  with  an  air  gap  in 
between.  The  walls  are  riddled  so  that  air  can  flow  in  the 
direction  opposite  to  the  main  heat  flow,  i.e.  from  outside  to 
inside.  Simulation  of  the  temperature  and  flow  fields  then 
allows  one  to  calculate  the  permeability  and  the  thermal 
conductivity  matrix  of  the  3D  structure.  With  this  informa¬ 
tion,  the  real  structure  can  be  replaced  by  a  two-dimensional 
“porous  material”  with  the  same  characteristics. 

This  is  shown  in  Fig.  11,  where  the  integrated  air  pre¬ 
heating  and  dynamic  insulation  of  an  SOFC  system  has  been 
simulated  in  lateral  direction.  The  outer  thermal  insulation 
has  been  modelled  as  a  porous  material  with  effective 
transport  parameters  that  were  obtained  from  the  detailed 
3D  model  shown  in  Fig.  10.  Furthermore,  details  of  the  stack 
have  been  neglected  here  since  for  testing  the  performance 
of  the  insulation  it  is  sufficient  to  specify  the  heat  and  mass 
fluxes  in  the  after-burning  zone  around  the  stack.  As  a 
consequence,  in  such  a  highly  optimised  2D  model,  typical 
calculations  of  the  coupled  pressure,  velocities,  and  tem¬ 
perature  fields  only  take  few  minutes  so  that  optimisation 
tasks  can  be  readily  performed. 


Fig.  9.  Performance  map  for  two  cells,  differing  in  the  electrolyte 
dimension:  thick  (top  graph)  and  thin  (bottom  graph).  Lines  of  constant 
efficiency  (solid),  electrical  power  (dashed)  and  fuel  utilisation  (dashed- 
dotted)  are  shown.  The  figures  for  the  current  are  given  in  A,  while  the 
others  are  dimensionless. 


Results  for  the  temperature  distribution  for  the  area 
enclosed  by  the  dashed  lines  are  presented  in  Fig.  12.  While 
there  are  large  temperature  gradients  within  the  heat 
exchange  zone,  the  temperature  becomes  rather  uniformly 
distributed  further  away  from  the  stack. 


air  gap 


heat  transport 


ceramic 

insulation 


Fig.  10.  The  3D  model  for  the  calculation  of  the  transport  properties  of  the 
thermal  insulation  structure. 
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Fig.  11.  The  2D  cross-sectional  model  of  a  SOFC  stack  with  integrated  air 
pre-heating  and  thermal  stack  insulation.  The  cold  air  enters  from  the  outer 
region,  passing  through  the  3D  structure  depicted  in  Fig.  10  and  enters  the 
heat  exchanger.  The  exhaust  gas  passes  through  the  heat  exchanger  and 
exits  through  the  outlets.  The  dashed  triangle  was  modelled  in  2D  using  a 
macro  transport  equation  for  the  thermal  insulation  region. 


Fig.  12.  Temperature  distribution  within  the  thermal  insulation  around  the 
SOFC  stack,  calculated  by  a  highly  efficient  2D  model.  The  shown  area 
corresponds  with  the  one  enclosed  by  the  dashed  lines  in  Fig.  11. 

4.  Conclusions 

4.1.  Summary 

We  have  presented  a  novel  approach  to  the  numerical 
modelling  of  fuel  cells,  stacks,  and  systems.  The  main 
feature  of  our  approach  is  the  adoption  of  macro  transport 
equations  to  describe  the  relevant  processes  in  terms  of 
averaged  physical  quantities.  The  method  of  numerical 
volume  averaging  exhibits  strong  analogies  to  the  work  of 
experimentalists,  in  the  sense  that  the  determination  of 
effective  parameters  can  be  described  as  performing  virtual 
experiments. 

It  is  possible  in  this  way  to  abstract  geometrical  details  of 
the  system  at  hand  and  to  absorb  their  influence  into 
effective  material  properties.  These  parameters  are  often 
tensor  quantities  reflecting,  e.g.  a  direction  dependent  heat 


conductivity.  This  enables  an  engineer  to  discretise  larger 
portions  of  the  device  at  hand  because  less  mesh  points  have 
to  be  used  in  an  averaged  domain. 

It  was  shown  that  a  hierarchy  of  models  can  be  set  up  for 
the  simulation  of  fuel  cells.  The  models  of  lower  levels 
provide  the  effective  parameters  for  the  next  higher  level.  It 
is  furthermore  possible  to  integrate  experimental  data  at  any 
stage  into  the  model  cascade,  which  either  is  used  to  validate 
the  model  or  to  input  otherwise  unknown  model  parameters. 
While  the  presented  examples  were  taken  from  SOFC 
systems,  the  method  applies  to  other  types  of  fuel  cells, 
such  as  PEFC  as  well. 

4.2.  Outlook 

There  is  a  lot  more  work  to  be  done.  On  one  hand,  it  is 
important  to  test  the  model  hierarchies  set-up  with  the 
NVAM  for  all  interactions  with  independently  measured 
experimental  data.  On  the  theoretical  side,  it  is  necessary  to 
find  rigorous  proofs  or  estimates  for  the  approximations 
made. 

Furthermore,  we  are  currently  applying  the  NVAM  to 
PEFC  systems. 
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